fOOL 


NPS55Gv75061 


/) 


NAVAL  POSTGRADUATE  SCHOOL 

Monterey,  California 


GAUSSIAN  APPROXIMATIONS  TO  SERVICE  PROBLEMS 

A  COMMUNICATION  SYSTEM  EXAMPLE 

by 

D.  P.  Gaver 
and 
J.  P.  Lehoczky 

June  1975 


Approved  for  public  release;  distribution  unlimited 
Prepared  for: 


FEDDOCS 
D  208.14/2: 
NPS-55GV75061 


NAVAL  POSTGRADUATE  SCHOOL 
Monterey,  California 

Rear  Admiral  Isham  Linder  Jack  R.  Borsting 

Superintendent  Provost 

The  work  reported  herein  was  supported  in  part  by  the  National 
Science  Foundation,  Grant  AG46  7,  at  the  Naval  Postgraduate  School, 
and  the  Air  Force  Office  of  Scientific  Research,  Grant  AFOSR74-2642 , 
at  Carnegie  Mellon  University. 

Reproduction  of  all  or  part  of  this  report  is  authorized. 

Prepared  by : 


UNCLASSIFIED 


SECURITY  CLASSIFICATION  OF  THIS  PAGE  (When  Data  Entered) 


REPORT  DOCUMENTATION  PAGE 


READ  INSTRUCTIONS 
BEFORE  COMPLETING  FORM 


1.     REPORT   NUMBER 

NPS55GV75061 


2.  GOVT  ACCESSION  NO 


3       RECIPIENT'S  CATALOG   NUMBER 


4.     TITLE  (and  Subtitle) 


Gaussian  Approximations  to  Service 
Problems:   A  Communication  System 
Example 


5       TYPE   OF    REPORT   A   PERIOD   COVERED 

Technical    Report 


6.     PERFORMING  ORG.   REPORT  NUMBER 


7.     AUTHORfaj 

Donald  P.  Gaver 
John  P.  Lehoczky 


8.  CONTRACT  OR  GRANT  NUMBER)-*.) 


9.  PERFORMING  ORGANIZATION  NAME  AND  ADDRESS 

Naval  Postgraduate  School 
Monterey,  California    93940 


10.     PROGRAM  ELEMENT,  PROJECT,   TASK 
AREA  ft  WORK  UNIT  NUMBERS 


11.     CONTROLLING  OFFICE  NAME   AND   ADDRESS 


12.     REPORT   DATE 

June    1975 


13.     NUMBER  OF   PAGES 


14.     MONITORING  AGENCY  NAME  4   ADDRESSf//  different  from  Controlling  Office) 


15.     SECURITY  CLASS,  (of  thie  report) 


Unclassified 


15a.     DECLASSIFI  CATION/ DOWN  GRADING 
SCHEDULE 


16.     DISTRIBUTION   ST  ATEMEN  T  (of  this  Report) 


Approved  for  public  release;  distribution  unlimited. 


17.     DISTRIBUTION  STATEMENT  (of  the  abatract  entered  In  Block  20,  If  different  from  Report) 


18.     SUPPLEMENTARY  NOTES 


19.     KEY  WORDS  (Continue  on  reverse  aide  It  necaaaary  and  Identify  by  block  number) 


Queues 

Service  Systems 


Diffusion 
Retrials 


20.     ABSTRACT  (Continue  on  reverae  aide  It  nacaaaary  and  Identity  by  block  number) 

Messages  arrive  at  a  group  of  service  channels  in  accordance 
with  a  time-dependent  Poisson  process.   An  arrival  either 
(i)  immediately  begins   k-stage  Markovian  service  if  an  empty 
channel  is  reached,  or  (ii)  balks  and  enters  a  retrial  popula- 
tion if  the  channel  sought  is  busy.   Diffusion  approximations 
to  the  number  of  messages  in  service  (each  stage)  and  in  the 
retrial  population  are  derived  by  writing  stochastic  differentia 


dd  ,; 


FaT73  1473 


EDITION  OF   1  NOV  65  IS  OBSOLETE 
S/N    0102-014-6601   | 


UNCLASSIFIED 


SECURITY  CLASSIFICATION  OF  THIS  PAGE  (When  Data  Entered) 


UNCLASSIFIED 


■,1_(-UW1TY  CLASSIFICATION  OF  THIS  PAGEfHTjen  Data  Entered) 


(1+0)   equations.   Steady-state  distributions  are  found  and 
compared  with  certain  simulation  results. 


UNCLASSIFIED 


SECURITY  CLASSIFICATION  OF  THIS  PAGEfWhsn  Data  Entered) 
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** 

J.  P.  Lehoczky 

Carnegie-Mellon  University 
Pittsburgh,  Pennsylvania 


1 .   Introduction. 

Any  reasonable  model  of  a  complex  communication  or  other 
service  system  requires  consideration  of  several  interacting 
populations.   Consequently,  it  becomes  necessary  to  name  the 
state  of  each  population,  and  to  describe  the  state  of  the  entire 
system  in  terms  of  a  vector-valued  state  variable.   With  rare 
exceptions  very  little  information  can  then  be  derived  in  a 
direct  manner,  e.g.  by  postulating  that  the  modelling  process  is 
multidimensional  birth-death  or,  more  generally,  Markov,  and 
deriving  mathematical  expressions  for  steady-state  probabilities, 
etc.,  as  exemplified  by  Feller  [3]  I,  Chap  XVII,  and  in  many 
papers.   Exceptions  do  occur  but  seem  rare,  see  the  cyclic  queue 
model  of  Gordon  and  Newell  [5],  and  related  work  by  Whittle  [10] 
and  by  Kingman  [6].   Consequently,  it  is  tempting  to  devise 
approximate  diffusion  models  for  such  processes;  previous  work  with 
this  intent  has  been  done  by  Schach  and  McNeill  [8] ,  McNeill  [7] , 
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and  the  approximation  studied  by  Barbour  [2].   Our  paper,  [4], 
indicates  how  this  may  be  done  for  a  transitory  situation.   In 
this  paper  we  study  a  communication  system  that  consists  of   c 
channels,  sought  for  by  arriving  messages  that  balk  if  they 
encounter  an  occupied  channel.   Balking  messages  enter  a  distinct 
(retrial)  population,  or  populations,  from  which  they  may  eventually 
either  defect,  or  again  apply  for  service  on  one  of  the  channels. 
Previous  work  in  this  problem  area  along  classical  point-process 
model  lines  is  described  by  Riordan  [9] .   Our  purpose  here  is  to 
indicate  how  stochastic  differential  equations  may  be  used  to 
write  down  a  model,  and  how  useful  information  may  then  be  derived. 
In  order  to  illustrate  the  degree  of  accuracy  achieved  we  have 
conducted  some  sample  simulations,  the  results  of  which  generally 
support  the  approximation  method. 


2.   The  Model. 

The  model  structure  is  depicted  in  Figure  1.   Input  to  a 
service  center  arrives  according  to  a  Poisson  process  with 
intensity   cA (t) .   The  service  center  consists  of   k   distinct 
compartments  and  a  total  of   c   channels  or  servers.   The  service 
process  on  the   i —  compartment  is  Markovian   (y  . (t) ) ;   i.e.  if 
a  "customer,"  here  message,  is  undergoing  stage-i  service  at   t 
it  terminates  within   (t,t+dt)   with  probability  y.(t)dt  +  o(dt), 
independently  of  previous  process  history.   Upon  completion  of  the 
i —  service  function  a  message  proceeds  immediately  to  the   (i+1) — . 
The  channels  are  considered  separate,  and  if  any  compartment  of  a 
channel  is  occupied  then  no  other  message  can  enter  that  channel. 
Consequently,  the  service  center  has  a  capacity  of   c   messages 
at  any  one  time,  and  a  single  message  has  a  total  service  time 
which  is  the  sum  of   k   independent  but  time-dependent  exponentials. 
If   y.(t)  =  y,   i=l,...,k   then  the  message  service  time  is 
Gamma   (k,y) .   In  fact,  the  compartments  are  introduced  in  order 
to  permit  the  modelling  of  non-exponential  service. 

We  postulate  that  a  message  selects  a  channel  uniformly 
at  random  from  among  the   c   possible  channels.   If  that  channel 
is  busy,  the  message  is  denied  immediate  access  to  service.   Other- 
wise, the  message  occupies  the  selected  channel  until  service  in 
all   k   compartments  has  been  completed.   The  random  selection  of 
channels  and  temporary  denial  of  service  goes  unrecognized  in  the 
context  of  classical  queueing  problems;  however, such  assumptions 
are  quite  appropriate  for  certain  communication  situations  (telephone 
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or  air-traffic  control,  for  example).   In  such  cases  the  customer 
may  not  automatically  select  an  empty  channel,  but  must  use  what- 
ever channel  is  appropriate,  whether  busy  or  not.   If  channels 
tend  to  be  used  equally,  the  random  selection  is  then  appropriate. 

If  a  customer  is  not  serviced  there  is  no  physical  queue. 
A  telephone  caller  receiving  a  busy  signal  can  only  hang  up,  per- 
haps with  the  intent  to  call  again  shortly.   An  airline  pilot  or 
submarine  attempting  to  contact  a  busy  traffic  controller  or  shore 
communication  facility  must  try  at  a  later  time,  possibly  to  be 
denied  service  again.   We  assume  that  messages  or  customers  who 
are  denied  service  enter  a  category  R.      Such  customers  may  retry 
or  recall  at  a  later  time.   These  customers  reapply  for  service 
after  being  in  state  R      for  a  Markovian   (v(t))~   period  of  time 
In  most  cases   v(t)   would  be  rather  large  if  the  customer  must 
receive  service  (air-traffic  control) ,  but  might  be  small  in  the 
case  of  telephone  calls. 

In  many  situations  customer  discouragement  will  occur, 

expecially  if  repeated  requests  for  service  are  denied.   We  adopt 

a  simple  model  of  customer  discouragement  by  assuming  a  binomial 

probability   1  -  a   that  a  customer  upon  leaving  R      decides  to 

leave  the  system  and  is  lost.   It  is  clear  that  a  vast  array  of 

models  can  be  introduced  to  represent  customer  discouragement. 

The  most  general  sort  might  include  states   iR. }     with   R. 

1  i=l 
signifying  that  the  customer  has  been  denied  service   i   times. 

The  holding  time  in  R .   is  Markovian   (v.  (t) )   ,   and   1  -  a .   is 

the  probability  of  loss.   We  discuss  only  the  simple  model 


mentioned  earlier  with  one  state   R,   but  the  reader  may  notice 
that  the  analysis  used  can  easily  be  extended  to  handle  the  most 

oo 

general  case  with  states   {R.  } 

We  wish  to  describe  the  behavior  of  this  system  by  calcu- 
lating system  characteristics  such  as  the  fraction  of  customers 
lost,  the  utilization  of  the   c   channels,  and  the  number  of 
customers  in  state   R,   all  as  a  function  of  the  system  parameters, 
Standard  Markov  chain  methods  can  be  used  on  this  model  owing  to 
the  Markovian  assumptions.   Unfortunately,  the  state  of  the  system 
is  a   k+2   dimensional  vector.   While  a  Markov  chain  analysis 
is  straightforward,  it  must  be  carried  out  numerically.   This 
makes  it  difficult  to  assess  the  influence  of  system  parameters 
on  the  quantities  of  interest.   Furthermore,  it  is  difficult  to 
deduce  information  about  the  transient  behavior  of  the  system  and 
to  handle  non-stationary  transition  probabilities  using  Markov 
chain  techniques.   For  these  reasons  we  adopt  the  method  of  diffu- 
sion approximations  as  an  approach.   This  consists  of  letting 
c  -y  +00   ancj  treating  the  resulting  random  processes  as  the  sum  of 
a  deterministic  process  plus  an  additive  noise  (diffusion)  process. 
This  technique  allows  one  to  compute  the  quantities  of  interest 
as  a  function  of  system  parameters  and  to  describe  the  transient 
behavior  of  the  system.   If  we  assume  the  transitions  are  station- 
ary, then  techniques  from  the  theory  of  stationary  processes, 
especially  spectral  analysis,  can  also  be  applied.   While  all 
results  are  exact  only  in  the  limit  as   c  -*■  °°,   we  present  compari- 
sons between  simulated  systems  with   c  =  10  and  20,  and  diffusion 


approximations.   The  comparisons  will  show  that  diffusion 
approximations  may  offer  surprising  accuracy  even  for   c   as 
small  as  10.   Thus  the  method  described  is  quite  appealing  as  an 
approximation  tool  in  the  present  problems ,  and  in  many  others 
as  well. 


3 .   Diffusion  Approximation  of  the  System. 
We  introduce  the  following  notation: 

Q. (t)  =  number  of  messages  or  customers  at  service  stage   i, 

at  time   t,   i  =  l,...,k, 

k 

Q(t)  =   £   Q.(t)  =  number  of  customers  receiving  service  at 

i=l   1 

time   t , 

R(t)  =  number  of  customers  in  state   R   at  time   t, 

L(t)  =  cumulative  number  of  customers  lost  by  time   t. 

A  customer  arriving  at  the  service  center  at  time   t   selects 
a  channel  at  random,  hence  with  probability  Q(t)/c   is  denied 
immediate  service,  and  with  probability   l-Q(t)/c   begins  service. 
It  is  straightforward  to  describe  this   k  +  2   dimensional  system 
and  its  transition  probabilities  over  the  interval   (t,t+dt); 
however,  we  choose  to  do  this  approximately  to  terms  of  order   dt 
using  techniques  from  the  theory  of  stochastic  differential  equa- 
tions; see  Arnold  [1],  and  also  Gaver,  Lehoczky  and  Perlas  [4]. 
Using  the  notation   dX(t)  =  X(t+dt)  -X(t)   for  a  stochastic  process 
{X(t),  t^O}   we  express  the  evolution  of  our  process  as  follows: 

dQ1(t)  =  (\(t)c+av(t)R(t))  (1  -Q(t)/c)dt- y1(t)Q1(t)dt 


+  /X(t)c(l  -  Q(t)/c)  dWA(t)  +  /av(t)R(t)  (1  -  Q(t)/c)  dWR(t) 

-  /y1(t)Q1(t)  dWQ  (t)  (3.1) 


dQi(t)  =  yi_1(t)Qi_1(t)dt-  yi(t)Qi(t)dt+  /yi_1(t)Qi_1(t)  dW     (t) 

i-1 

-  /y . (t)Q±(t)  dW   (t)    for   i  =  2,...,k 
dR(t)  =  -(l-a+a(l-Q(t)/c) )v(t)R(t)dt + X(t)c(Q(t)/c)dt 


-  /(l-a+a(l-Q(t)/c))v(t)R(t)  dWR(t)  +  A  (t)c  (Q(t)/c)  dWA  (t) 


dL(t)  =  (l-a)v(t)R(t)dt+  /(l-a)v(t)R(t)  dWR(t)  . 


In  equations  (3.1)  the  terms   dW   (t) ,   dW  (t)  ,   and 

dW  (t)   are  the  "derivatives"  of  independent  standard  Wiener 
A 

processes,  and  as  such  are  usually  called  Gaussian  white  noise. 
Other  approximations  are  possible,  but  are  not  pursued  here.   We 
mention,  for  example,  using  Poisson  white  noise  where  the  standard 
Wiener  process  is  replaced  by  a  Poisson  process  with  zero  drift. 
A  word  about  the  derivation  of  our  equations  (3.1)  is  in 
order.   Consider  that  for  the  occupancy  of  the  first  service 
compartment,   Q, (t) :   conditional  on  the  values  of   Q, (t)   and 
R(t)   the  drift  or  mean  change  of   0,   in  time   dt   is  (a)  positive 
by  the  amount   (A  (t)  c+ctv  (t)  R(t)  )  (1  -     )dt,    where  the  first 
term  represents  the  expected  number  of  arrivals  in   (t,t+dt), 
and  the  other  represents  the  probability  of  acceptance  into  ser- 
vice in  compartment   i =  1 ,   while  (b)  the  second,  negative,  term 
-]ii  (t)Q,  (t)dt,   represents  the  expected  number  departing   Q,   in 
(t,t+dt);   hence  the  difference  is  the  net  expected  increase  in 
Q, .   Now  for   dt   small  we  represent  the  fluctuating  (diffusion) 


component  of  input  by  (c) :  /A  (t) c(l  -  Q (t)/c)  dW^(t),   where  the 


scale  factor  is  the  standard  deviation  of  a  Poisson  process  with 
mean   X (t)c (1  -  Q (t)/c)   and   dW^ (t)   is   W(0,dt),   plus  (d)  a 
corresponding  but  independent  term  representing  arrivals  from   R, 
minus  (e)  another  corresponding  term  representing  departure  from 
Q, .   Although  this  derivation  is  heuristic,  the  rationale  is 
simply  that  in  an  interval  of  length  dt   the  various  condition- 
ally Poisson  components  of  change  are  approximately  normal, 
owing  to  the  fact  that   c   is  presumed  large. 

We  wish  next  to  view  the  state  of  the  system 
(Q-.  (t)  ,  .  .  .  ,Q,  (t)  ,R(t)  ,L  (t)  )   as  the  sum  of  a  deterministic  pro- 
cess plus  a  noise  or  diffusion  process.   To  accomplish  this  we 
introduce  the  standardized  noise  variables: 

Q.  (t)  -cq.  (t)  k 

X,  (t)  =  -± i  =  l,...,k;   X(t)  =  I      X.  (t) 

/c  i=l 


Y(t)  =  R(t)  -  cr(t) 


z(t)  =  L(t)  -cMt) 


-  c 
/c" 


(3.2) 


k 
q(t)  =  I      q  (t) 
i=l 

We  then  write   (Q, (t) , . . . ,Q,  (t) ,R(t) ,L (t) )   in  vector 


fashion  as 

c(q1(t)  ,q2(t)  ,  ...,qk(t)  ,r(t)  ,  fL  (t)  ) 

+   c(X1(t) ,X2(t) ,... ,Xk(t) ,Y(t) ,Z(t) ) , 
the  first  term  is  the  deterministic  approximation,  the  last  is 


10 


the  noise  process.   We  substitute  definitions  (3.2)  into  (3.1) 
and  divide  the  resulting  equations  by   /c.   The  results  are,  to 
terms  of  order  1  in   c,   expressible  as 

dX1(t)  =  {(1  -q(t))av(t)Y(t)  -  (X(t)+av(t)r(t)  )X(t)  -  y1(t)X1(t)  }dt 


+  /X(t)  (1  -q(t))  dWx  +  /av(t)r(t)  (1  -q(t)  )  dWR  -  /y1(t)q1(t)  dWQ 


•c{q{(t)  -  (X(t)+av(t)r(t))  (1  -q(t))  +y1(t)q1(t)  }dt, 


dxi(t)  =  {yi_1(t)xi_1(t)  -  yi(t)xi(t)}dt  +  /yi_1(t)qi_1(t)  dW 

i-1 


-  v^ITTtTqTTty  dW   -  /c"{q[(t)-yi_1(t)qi_1(t)  +  yi(t)qi(t)  }dt 

i 

for   i  =  2,...,k,    (3.3) 
dY(t)  =  {-[(l-aq(t))v(t)]Y(t)  +  [X(t)  +av(t)r(t)]X(t)}dt 


-  /(l-aq(t))v(t)r(t)  dWR  +  /X(t)q(t)  dW^ 


-  /c"{r'  (t)  +  (l-aq(t))v(t)r(t)  -  X(t)q(t)  }dt, 
and 

dZ(t)  =  (l-a)v(t)Y(t)dt +  /(l-a)v(t)r(t)  dWR 

-  /c"{£'(t)  -  (l-a)v(t)r(t)  }dt. 

We  now  let   c  -»■  °°   in  equations  (3.3).   Clearly,  in  all 
cases  the   /c   term  must  be  identically   0,   or  else  the  equations 

11 


explode.   Setting  the  Jc      term  to   0   we  derive  a  system  of 
ordinary  differential  equations  satisfied  by  the  deterministic 
approximation : 

q|(t)  =   (X(t)+av(t)r(t))  (l-q(t))  -  y1(t)q1(t) 

q|(t)  =  Ui_1(t)qi_1(t)  -  u±(t)qi(t)      i=2,...,k       (3.4) 

r'(t)  =  -(l-aq(t))v(t)r(t)  +  X(t)q(t) 

V  (t)  =  (l-a)v(t)r(t) 

It  is  easy  to  find  q.  (t)   in  terms  of  <3-_i  (t)   for 
i  =  2/.../k   by  solving  the  second  equation  in  (3.4). 


q.tt) 


ft   -/sy.(x)dx  -/0vi.(s)ds 

e  yi_1(s)qi_1(s)ds  +  qi(0)e 


0 


but  no  further  explicit  results  seem  attainable  without  simplifi- 
cation and  further  parameter  specification.   Of  course  numerical 
solution  of  the  differential  equations  by  computer  is  always 
possible,  and  may  well  prove  to  be  the  fastest  route  to  useful 
information. 

Steady  State  Behavior  of  the  System 

Suppose  however  that  we  let   t  -*■  °°  and  attempt  to  find  a 
steady  state  solution.   As   t  -*  °°   let   X  (t)  ,   v(t),   y.(t), 
q.(t),   and   r(t)   converge  to   X,   v,  ]i  •  r      q-/   and   r   respec- 
tively.  Then  q! (t)   and   r*  (t)   all  converge  to   0   and  the 
steady  state  equations  become 

12 


0  =  (X+avr)  (1-q)  -  y.^ 

0  =  yi-lqi-l  "  yiqi      i  =  2"-"k 
0  =  -(1-aq) vr  +  Xq 

£•  (oo)  =  (l-a)vr 

k  , 
Letting   y  =  1/  £  —  we  find  the  steady-state  values 

i=lyi 

r(co)  =  Xc3  A+y-/(X+y)^-4Xay 

*  ;    v(l-aq)     q^  ;  2ay 

The  negative  square  root  is  used  for   q   since   0  £  q  £  1.   Fur- 
thermore/ the  asymptotic  loss  rate,  I'  (<*>)  ,         is  given  by 

— — -.   The  input  rate  is   X,   so  the  output  rate  of  those 

1-aq  r 

actually  served  is   X  (..  J"°  )  ,   and  the  asymptotic  loss  fraction 
(the  fraction  of  customers  leaving  without  receiving  service)  is 
given  by   (1-a) q/ (1-aq) . 

We  intend  to  carry  out  the  subsequent  analysis  under  the 
assumption  that  steady  state  conditions  prevail.   In  this  model 
the  deterministic  equations  are  nonlinear,  hence  there  is  a  ques- 
tion of  the  stability  of  the  system  (in  the  sense  of  Liapunov) . 
We  must  be  concerned  with  the  effect  of  a  small  perturbation  when 
the  system  is  in  steady  state.   If  the  system  is  not  stable  then 
it  will  diverge  from,  rather  than  return  to,  steady  state.   The 
stability  can  be  established  when   X(t),   y. (t)   and   v(t)   con- 
verge to   X,   y.,   and   v   by  first  linearizing  equations  (3.4). 
We  omit  l(t)       from  the  system  since  it  clearly  grows  without 


13 


(3.5) 


limit  and  has  no  steady  state  value.   We  assume  now  that  the 
parameters   X,   y.,   and   v   are  all  constants,  and  express  (3.4) 
in  matrix  form: 


q{(t) 
qj(t) 


qi<t) 

r'(t) 


r-(X+u1) 


■X    av-\  /-q1  (t)  -\ 


0 


k 
X    -v 


qk(t) 
r(t) 


(3.7) 


or 


X1 (t)  =  c  +  AX(t) 


The  system  will  be  stable  in  the  sense  of  Liapunov  if  the 
k  +  1   eigenvalues  of  the   A  matrix  have  strictly  negative  real 
parts.   The  eigenvalues  all  have  strictly  negative  real  parts 
provided  y.  >  0  and  v  >  0   for   i  =  l,...,k.   We  prove  this  in 
the  appendix. 

We  now  turn  to  a  description  of  the  noise  process.   We 
assume  deterministic  equations  (3.7)  are  satisfied.   The  represent 
tation  (3.3)  becomes  in  matrix  form 


where 


dU(t)  =  A.U(t)dt  +  B.dW. 

r~  ^t'*'  *^t  **'t 

U(t)  =  (Xx(t) ,... ,Xk(t) ,Y(t) ,Z(t)) ' 


W(t)  =  (Wx(t)fWR(t) ,WQ  (t),...,WQ  (t))' 


(3.8) 
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K  = 


f- (y, +X+avr) 


- (X+avr) 

y2 
0 


0 

0 

X+avr 

X+avr 

0 

0 

-(X+avr)    av  (rq)    0^ 
0  0       0 


0 

-Vk 

0 

X+avr 

- (1-aq) v 

0 

(l-a)v   0 

and 


st  ■ 


r/X (1-q)    /avr (1-q) 
0  0 


0 

0 

/Xq 

-/(1-aq) vr 

0 

/(1-aq) vr 

/yq 

/yq 

0 


0 

/yq 

/yq" 

0 


0  > 


0 

•/yq 

0 
0 


In  the  above  definitions  of   A^   and   B,_   we  have  omitted  the 

~t       ~t 

argument   t   in   X  (t)  ,      y.(t),   v(t),   q(t),   and   r(t)   for 

notational  convenience.   As   t  •*■   <»  a.   and   B.   converge  to   A 

'^t       ~t  ~ 

and   B   again  given  by  (3.8);  however,  asymptotically   Z(t)   is  not 
of  interest  because   L(t)  ■*■   +».   Furthermore,   A,,   is  singular. 
For  these  reasons  we  eliminate   Z(t)   and  reduce  the  dimension  to 
k  +  1.   The  resulting  stochastic  differential  equation  becomes 
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dX<t>  =  C^vmdt  +  J^d^  (3.9) 

with  V(t)  =  (X-,  (t)  , ..  .,Xk(t)  ,Y(t))  '   and   C   given  by  Afc  with 

the  last  row  and  column  removed.   D,   is  given  by   B.   with  the 

~t  J      ~t 

last  row  removed.   Again   C.   and   D,   converge  to   C   and   D   as 

-'     *^t       *^t  ***  **" 

t  ■*■   °°. 

The  stochastic  process  V(t)   satisfying  (3.9)  is  a  non- 
stationary  multivariate  Ornstein-Uhlenbeck  process.   This  process 
has  been  extensively  studied,  with  many  results  recorded  by  Arnold 
([1],  Section  8.2).   We  shall  make  use  of  these  results  in  what 
follows. 

We  may  integrate  (3.9)  to  find 


V(t)  =  V(0)  +    C4_V(t)dt  + 


t 


t 
D.  dW^..  (3.10) 

o  ^  T- 


The  process  V(t)   is  Gaussian  if   V(0)   is  either  constant 
or  itself  Gaussian.   This  is  clear  from  direct  examination  of 
(3.10).   Suppose   V(t)   is  either  constant  or  Gaussian.   V(t+dt) 
is  the  sum  of  V(t)  +  C.V(t)dt,   which  is  Gaussian,  and  another 
Gaussian  variable.   Thus   J/(t+dt)   will  also  be  Gaussian  and  it 
remains  to  characterize  the  marginal  moments  of  V(t)   as  well  as 
the  covariance  structure. 

Let   u^.  =  E(V(t)),  lt  =   E((V(t)-jj#t)  (V(t)-Jit))  '    be  the 
marginal  mean  and  covariance  of  V(t).   It  is  shown  in  Arnold  [1] 
that  u,   and   £.   satisfy  first-order  differential  equations. 
First,  the  mean  vector  is  described  by 

it  =  £t*t     with     k*  =  ^(0)  (3-11} 
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and,  second,  the  covariance  matrix   £    is  the  unique  symmetric 
nonnegative  definite  solution  of  the  matrix  differential  equation 

It  =  £tlt  +  It£i  +  fit&    with    lo  =  E(^(0)-Jio)(^(0>-iio),)-   (3-12) 

Equations  (3.9),  (3.11),  and  (3.12)  can  be  formally  solved 
with  the  aid  of  the  fundamental  matrix   $(t),   that  is  the  matrix 
of  solutions  of  the  homogeneous  equation   $(t)  =  C,$(t),   $(0)  =  I. 
For  example,  if   C   =  C,   in  steady  state,   $(t)  =  exp(Ct).   Using 
£(t)   we  find 

ft 


V(t)  =  $(t)  (V(0)  + 


($(s))"1D  dW  ) 
0  ~      ~s  ~s 


Et  =  £(t)E(V(0)) 

(3.13) 

K(s,t)  =  EdVis)-}^)   (V(t)-U^)  ') 

min (s, t) 

o 

letting   s  =  t,   k(s,t)  =  £    thus  providing  a  formula  for   £. . 
We  note  that  in  practice  it  will  often  be  convenient  to 
apply  a  computer  routine  for  solving  first-order  differential 
equations  directly  to  (3.11)  and  (3.12). 

Suppose  now  we  assume  C,_   =   C      and   D,_  =  D,   that  is  we 

~t   ***       *wt   **• 

are  in  steady  state.   It  will  be  shown  in  the  appendix  that  all 
k  +  1   eigenvalues  of   C   have  negative  real  parts.   In  fact,  the 

matrix   C   is  nearly  identical  to  the  matrix  A  examined  earlier. 

'°°   Ct      C't 
Then  if   v(0)  ~N(0,7)   with   7  =    e~  DD'e~  dt ,   V(t)   is  a 

Jo 

stationary  Gaussian  process  with   E(V(t))  =  0,   E (V (t) (V (t) ) ' )  =  7 
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where   £   is  as  defined  above  or  is,  equally,  the  unique  nonnegative 
definite  solution  of  the  equation 

£i  +  l£'  =  -dd' . 

Furthermore,  the  covariance  function   K(s,t)  =  H(s-t)   is  given 
by 


C(s-t) 


£    s  £  t  £  0 


(3.14) 


£e£,(t~s)    t  *s  *  0  . 


We  summarize  our  description  of  the   k  + 1   dimensional 
queueing  system  as  follows  using  the  diffusion  approximation: 

(Q1(t)  /...,Q]c(t)  ,R(t))  *r 

C(q1(t) ,... ,qk(t) ,r(t))  +  /C (X1 (t) ,. .. ,Xk(t) ,Y(t)) 

where  the  first  term  is  given  by  (3.4)  and  the  second  is  a  multi- 
variate Ornstein-Uhlenbeck  process  with  mean  jj^  ,   covariance 
function  H(s-t)   as  described  earlier. 

Results  for  the  Single  Service  Compartment 

We  give  the  exact  formulas  in  the  steady  state  case  for 
the  special  case  of   k  =  1,   a  single  service  compartment.   In 
this  case,  the  deterministic  approximations  are  still  given  by 
(3.6).  The  noise  approximation  will  be  a  bivariate  Ornstein-Uhlen- 
beck process  with  mean   0.   The  covariance  matrix   7   will  be  the 
unique  symmetric  nonnegative  definite  solution  of  the  equation 
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Al    +  lA'  =  -BB 


*••*"      t**r~ 


(3.15) 


where 


A  = 


B  = 


C-(y+X+avr)     (l-q)av 
X  +  avr     - (1-aq) v 
/X(l-q)   Axvr (1-q)   -/y 
/Xq    -/(l-aq)vr    0 


^ 


BB'  = 


bl   b3 


Lb_   b„ 


for  notational  simplicity  in  what  follows.   Lastly 


I  =  (°ij)      ail  =  Var(x(t)) 


a22  =  Var(Y(t)) 

a12  =  Cov(X(t) ,Y(t)) 


Solving    (3.15)    we    find 


a..    =    [b,  (  |A|+a;o)  -  2a10a00b0  +  a'bJ/D 


11 


l22 


12  22  2    12  3 


°12  =  [-a21a22bl  +  2b2alla22"b3alla12]/D 


a        =    [a*  b,  -  2a,  ,a9  b9  +  h  (  lAl+a,2,  ]/D 


22 


■21"1   ""11"21"2   "3 


11 


D  =  2(ai;L+a22)  |a|  . 


(3.16) 


In  the  last  section  we  present  a  comparison  of  the  theoretical 
calculations  and  simulation  results  in  a  variety  of  situations 
with   c  =  10  and  20.   We  are  approximating  the  state  of 
the  system   (Q(t),R(t))   by   c (q, r)  +  /c (X (t ) , Y (t) ) .   The  diffusion 
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approximation  allows  for  more  than  a  description  of  the  marginal 
behavior  in  that  the  transient  behavior  can  be  characterized  and 
the  joint  behavior  at  any  set  of  time  point,   (t, ,...,t  )   can 
be  worked  out  using  the  multivariate  Ornstein-Uhlenbeck  process. 
However,  we  do  not  present  numerical  results  of  this  kind  in 
this  paper. 
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4 .   The  Spectral  Matrix. 

We  now  assume  the  system  is  in  steady  state  and  adopt  the 
steady  state  diffusion  approximation  of  it.   The  noise  process  is 
given  by 

dU(t)  =  AU(t)dt +  BdW(t) .  (4.1) 

The  noise  process   {U(t),t^0}   is  stationary,  hence 
(Q-.  (t)  ,  .  .  .  ,R  (t)  )   will  also  be  stationary.   This  observation  opens 
up  the  possibility  of  applying  techniques  from  the  theory  of 
stationary  processes  to  our  queueing  process.   In  particular  we 
compute  the  spectral  matrix  associated  with  the  noise  process. 

We  begin  with  the  spectral  representations  of  the  stationary 
processes   U(t)   and  W(t) 


U(t)  = 


ela)t  dS0(co)  ,    W(t)  =    ela)t  dSw(oj)         (4.2) 


where   {STT(w),w  e  (-co,<»)  }   and   {Sn(o)),iDe(-»,<»)}   are  processes  with 

orthogonal  increments.   The  spectral  matrices  associated  with   U 

and  W,   f„   and   frT,   are  the  matrices  of  cross  spectral  densities 
~   ~u       ~W 

of  the  stationary  processes  and  are  given  by 
(w) 


£7   =  E(dSu(o))dSu(a>)),   fwU)  =  E(dSw(a))dSw(W))  =  ^. 

We  wish  to  compute   fTT(w)   in  terms  of   A  and   B. 

Combining  (4.1)  and  (4.2)  by  differentiation  of  (4.2)  we 
find 

(iwl-A)dSTT  =  BdS..  (4.3) 
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Taking  transposes  in  (4.4),  multiplying,  and  taking  expectations 
we  obtain 


BB' 


E  [  (iu)I-A)  dSyd^  (iu)I-Aj  ]  =  E  [Bd^dSwB  ']  =  ^ 


^*  *w 


rO  *^- 


Solving    (4.5)    for      frT   =   E(dSTTdSTT)      we    find 

fjfjM    =   ^(A-ioilj'^^A'+io)!)"1. 

^» 

We  illustrate  this  computation  in  the  case   k  =  1,   a 
single  service  compartment.   Here 


A  = 


-(y+X+avr)    (l-q)av^ 


,   BB' 


2yq    A  ^ 
A    2\q) 


K      X  +  avr     -(1-aq)  \)J 

with   A  =  X/q(l-q)  -  Vr/a(l-q)  (1-aq) . 

After  a  few  simple  matrix  inversions  and  multiplications  we 

compute 


£u("> 


rfQQ(w)   fQR(w)> 


LfRQ(w)   fRR(w)J 


(4.4) 


(4.5) 


where 


fQQU)  =  [2yqoa2  +v  (1-aq)  (2yqv  (1-aq) +Aav(l-q)  ) 


+  av(l-q)  (Av(l-aq)  +  2Xqav(l-q) ] /D 


fRR(to)  =  [2Xqw2 +  (X+avr) 2yq(X+avr+A (y+X+avr) ) 


+  (y+X+avr)  ( A(X+avr)  +  2Xq (y+X+avr) ] /D     (4 . 6) 
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fQR(u>)  =  {  [Aw2  +  v(l-aq)  (2yq(X+avr)  +  A  (y+X+avr)  ) 


+  av (1-q)  (A (X+avr)  +  2 Xq (y+X+avr) ) ] 
+  i(0  [AV  (1-aq)  +  2Xqav  (1-q)  -  2yq  (X+avr) 

-  A (y+X+avr) ] }/D 


fRQ(a))  =  fQR(a,) 


with 


D  =  [a)2  +  av  (1-q)  (X+avr)  -  X  (1-aq)  (y+X+avr)  ]  2 

+  [V  (1-aq)  +  y  +  X  +  avr]  2a>2  . 

The  functions   f0Q(u)),   f   (go)   are  the  spectral  densities 
of   Q(t)   and   R(t)   respectively.   The  real  part  of   f__(u)) 
is  the  cospectral  density  and  the  imaginary  part  is  the  quadrature 
spectral  density.   The  latter  two  densities  provide  information 

about  the  phase  behavior  of   (Q(t),R(t)).   Notice  that   f   (w) 

-2 
and   fR  (w)   exhibit  tail  behavior   w    and  thus  fluctuate  in  a 

high-frequency  manner  similar  to  that  of  the  ordinary  Ornstein- 

Uhlenbeck  process. 
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5.   Comparison  with  Simulation 

The  previously  quoted  results  are  derived  under  the  supposition 
that   c,   the  number  of  service  channels,  becomes  indefinitely 
large.   Since  it  will  be  desirable  to  apply  the  approximations  in 
case   c   is  finite,  we  have  undertaken  to  perform  several  simula- 
tions of  particular  systems  with  moderate   c-values,  and  thus  to 
provide  an  empirical  check  of  model  accuracy.   We  shall  see  that 
our  approximations  are  generally  good. 

The  simulation  simulated  is  the  system  (Q (t) ,R(t) ,L (t) )  in 
which  k  =  1   (the  single  service  compartment  case),  and   X,   u, 
and   v  are  constants.   By  virtue  of  the  Markov  nature  of  the 
system  we  see  that  (i)  sojourns  in  states  are  of  independent  and 
exponentially  distributed  duration,  and  (ii)  state  changes  occur 
in  accordance  with  multinomial  Bernoulli  trials.   Reference  to 
(  3.1  )  shows  that  the  state-dependent  transition  rates  are  the 
following 


(Q(t) ,R(t) ,L(t))  +  (Q(t)+l,R(t) ,L(t)) 

+  (Q(t)+l,R(t)-l,L(t) ) 

■*  (Q(t)-l,R(t)  ,L(t)+l) 

-*  (Q(t)  ,R(t)+l,L(t)) 

-v  (Q(t)  ,R(t)-l,L(t)+l) 


Xc[l-Q(t)/c] 

VCtR(t)  [l-Q(t)/c] 

uQ(t) 

XcQ(t)/c 

v(l-a)R(t) 


(5.1) 


Hence,  given  that  at  time   t   the  system  is  in  state   (Q (t) ,R(t) ,L (t) 
it  resides  there  for  independently  and  exponentially  distributed 
times  with  means  equal  to  the  inverse  of  the  sum  of  the  right-hand 
side  of  (5.1),  and  then  instantaneously  jumps  to  a  new,  neighboring, 


24 


state  with  probabilities  given  by  the  right  side  of  (5.1)  times 
the  mean  sojourn  time  in  state.   Our  simulation  program  is  based 
on  this  scheme. 

The  simulation  results  recorded  were  the  fractions  of  times 
that  the  system  inhabited  each  state   (Q=i,  R= j ,    L=k)   over  the 
course  of  the  simulation,  the  latter  chosen  to  be  long.   These 
fractions  approximate  the  stationary  state  probabilities.   The 
latter  probability  estimates  were  in  turn  used  to  compute  estimates 
of   E[Q(oo)],   Var[Q(°°)],   E[R(«>)],   Var[R(»)],   and   E[L(°°)].   Also, 
the  tabulated  empirical  marginal  distributions  of  simulated   Q 
and   R  were  plotted  on  normal  probability  paper  in  order  to  pro- 
vide a  visual  check  for  normality. 

Discussion 

Agreement  of  the  diffusion  model  with  the  simulation  is, 
apparently,  quite  acceptable  for  the  cases  considered,  which  by 
design  include  relatively  small  numbers  of  channels.   Extensive 
additional  simulation  results,  left  unreported  here,  convey  the 
same  message.   The  most  noticeable  discrepancy  occurs  in  the 
estimates  of   Var[R(°°)]:   that  obtained  from  the  analytical 
approximation  consistently  exceeds  the  simulation  estimate.   Attempts 
to  show  that  simulation's  failure  to  reach  steady  state  (starts 
were  normally  made  at   Q(0)  =R(0)  =  0)   by  starting  higher  gave 
essentially  the  same  results.   The  discrepancy  remains  unexplained. 

In  order  to  aid  quick  appraisals  we  have  tabulated  some  key 
quantiles;  on  the  simulation  side  these  are  inexact  both  because 
of  the  inherent  discreteness  of  the  distributions  and  because  of 
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simulation  sampling  error,  and  on  the  diffusion  side  because  of 
the  use  of  the  continuous  normal  approximation  to  a  discrete 
distribution.   We  have,  for  example,  diffusion  approximated 

Q0.95   bY   ECQ]  +  1#65  St-Dev-[Q1'   Qo.75   by   E[Q0.75]  + 
0.68  Std.Dev. [QQ  75]  ,   and  QQ  5Q   by  E[Q].   The  selected  proba- 
bilities were  calculated  using  a  simple  continuity  correction, 
e.g. 


P{Q  s:  x}  =  — 


"2tt 


(x+|-E[Q]) /Std.Dev. [Q] 

-  |z2dz 
e 


The  latter  give  at  a  glance  an  impression  of  the  Gaussian  marginal 
model  adequacy,  which  by  and  large  is  good  out  to  the  two-sigma 
level.   As  is  to  be  anticipated,  the  Gaussian  approximation 
degenerates  in  quality  when  E[Q]  +  k  Std.Dev. [Q]   nears  either 
zero  or   c. 
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Table  1 
Comparison  of  Simulation  and  Diffusion  Approximation 
X   =  7,   u  =  4,   v  =  6,  a   =    0.5,   c  =  10 


EQQ(oo)] 
Var[Q(°°)  ] 
Std.Dev.  [Q(°°)  ] 
E[R(»)] 
Var [R(») ] 
Std.Dev.  [R(°°)  ] 


Simulation 

7.31 

2.00 

1.41 
13.55 
16.84 

4.10 


Diffusion  Approximation 


7.34 
2.08 
1.44 
13.54 
20.80 
4.56 


Q 


Q 


r0.05 
0.25 
!0.50 
!0.75 
0.95 
'0.05 
'0.25 
'0.50 
'0.75 
'0.95 


4 

6 

7 

8 

9 

7 

10 

13 

16 

20 


5 

6 

7 

8 

10 

6 

10 

14 

17 

21 


P{Q(« 
P{Q(< 
P{Q(< 
P{R(< 
P{R(< 
P{R(' 


£9} 
£7} 
*  5} 
£  20} 
£  13) 
£.6} 


956 
531 
104 
946 
518 
030 


933 
544 
102 
937 
496 
062 
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Table  2 
Comparison  of  Simulation  and  Diffusion  Approximation 
X  =  5,   y  =  4,   v  =  6,  a   =    0.5,   c  =  10 


Simulation 

Diffusion  Approximation 

E[Q(»)] 

6.45 

6.49 

Var[Q(°°)] 

2.36 

2.51 

Std.Dev.  [Q(°°)  ] 

1.54 

1.58 

E[R(«)] 

8.09 

8.01 

Var[R(°°)  ] 

10.20 

12.65 

Std.Dev.  [R(°°)  ] 

3.19 

3.56 

Q0.05 

3 

4 

Q0.25 

5 

5 

Q0.50 

6 

6 

Q0.75 

7 

8 

Q0.95 

8 

9 

R0.05 

3 

2 

£\.  **.   «  .- 

5 

6 

0.25 

R0.50 

8 

8 

p 

10 

10 

0.  75 

R0.95 

13 

14 

p{Q(oo)  s:  9} 

.986 

.971 

P{Q(oo)  *  6} 

.497 

.504 

p{Q(oo)  <:  3} 

.030 

,030 

P{R(°»)  £  14} 

.966 

.965 

P{R(°°)  £  8) 

.581 

.555 

P{R(°°)  £  2} 

.021 

.050 
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Table  3 
Comparison  of  Simulation  and  Diffusion  Approximation 
A  =  7,   u  =  4,   v  =  6,  a   =    0.5,   c  =  20 


Simulation 

Diffusion  Approximation 

E[Q(«)] 

14.62 

14.69 

Var[Q(°°)] 

4.08 

4.16 

Std.Dev.  [Q(°°)  ] 

2.02 

2.04 

E[R(«)] 

27.01 

27.08 

Var  [R(°°)  ] 

34.07 

41.60 

Std.Dev. [R(~) ] 

5.84 

6.45 

Q0.05 

11 

11 

Q0.25 

13 

13 

Q0.50 

14 

15 

Q0.75 

16 

16 

Q0.95 

17 

18 

R0.05 

17 

16 

R0.25 

22 

23 

R0.50 

26 

27 

R0.75 

30 

31 

R0.95 

36 

38 

p{Q(oo)  i  19} 

.998 

.991 

p{Q(oo)  s:  15} 

.654 

.655 

P{Q(co)  &   11} 

.066 

.059 

P{R(°°)  ^  40} 

.982 

.981 

P{R(<»)  *  27} 

.532 

.526 

P{R(~)  ^  14} 

.026 

.025 
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Table  4 
Comparison  of  Simulation  and  Diffusion  Approximation 
X    =    5,   y  =  4,   v  =  6,  a   =  0.5,   c  =  20 


Simulation 

Diffi 

usion  Approximation 

E[Q(°°)] 

12.89 

12.98 

Var[Q(°°)  ] 

4.89 

5.02 

Std.Dev. 

[Q(»)l 

2.21 

2.24 

E[R(«)  ] 

15.98 

16.02 

Var  [R(°°)  ] 

20.58 

25.30 

Std.Dev. 

[R(« 

>)] 

4.54 

5.03 

Q0.05 

9 

9 

Q0.25 

11 

11 

Q0.50 

12 

13 

Q0.75 

14 

15 

Q0.95 

16 

17 

R0.05 

8 

8 

R0.25 

12 

13 

R0.50 

15 

16 

R0.75 

18 

19 

R0.95 

23 

24 

P{Q(«)  * 

17} 

.989 

.978 

P{Q(oo)  ^ 

13} 

.596 

.591 

P{Q(oo)  ^ 

9} 

.066 

.061 

P{R(«>)  ^ 

24} 

.960 

.954 

P{R(«>)  £ 

16} 

.569 

.540 

P{R(°°)  ^ 

8} 

.037 

.067 
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Appendix 

We  prove  that  the  matrices  £   from  (3.7)  and  £    from 
(3.9)  have  eigenvalues  with  strictly  negative  real  parts.   Both 
matrices  have  the  form 


M  = 


-(U-j+c) 

-c 

-c 

-c 

yi 

-y2 

0 

0 

0 

U2 

-v3       .. 

•                          • 
• 

• 

0 

y3 

• 

• 
• 

• 
• 

0 

0 

0 

* 

0 

0 

-V 

c 

c 

c 

c 

o 


cJ 


with   c  =  X,   d  =  av,   e  =  v   for  &  and   c  =  A+avr,   d  =  av(l-q)  , 
e  =  v(l-aq)   for   $L_.   In  both  cases   e-d  =  v(l-a)  >  0. 

We  wish  to  solve  for  the   k  +  1   roots  of  the  equation 
|M  -  61  I  =  0.   Simple  manipulation  gives 


M-  81   = 


(8+y^ 


1 
0 


0 
c 


-(9+y2) 


V 


2 
0 

* 

0 

c 


(6+yk) 


(8+e-d) 

0 


0 
-(8+e) 


k+1 


k+3 


=  (~1)K  X  (8+e)  7T  (8+y.)  +  (-1)"  J  (8+e-d)Dk 

i=l     1 
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where 


D,  =  c 


(6+y2) 


(9+Vl» 


"k-1 


-(6+yk) 


The  determinant  D,   can  be  computed  recursively  as 

Dk  ■  (6+VDk-i +  ■    "i- 

1=1 
and  these  equations  can  be  solved.   We  find 


Dk  = 


n  (Q+M±)    -  n  yi 
i=l         i=l 


/e 


(A.l) 


Substitution  of  (A.l)  into   I M  —  01 1   gives 


0  =  (8+e)  II  (6+y.)  +  c(G+e-d) 
i=l     1 


k  k    > 

n  (e+y.)  -  n  y.  /e. 

i=l     X    i-1  ^ 


(A. 2) 


We  wish  to  show  that  all   k  +  1   roots  of  (2)  have  strictly 
negative  real  parts  assuming  e,   e-d,   y,  ,  .  .  .  ,   y,   are  all 
strictly  positive. 

The   k  +  1   degree  polynomial  on  the  right  side  of  equation 
(A. 2)  has  strictly  positive  coefficients,  thus  by  Descartes  rule 
only  strictly  negative  real  roots  are  possible.   It  remains  to 
consider  the  case  of  complex  roots  and  we  assume   9  -  a  +  b..   Since 
complex  roots  of  (A. 2)  appear  in  conjugate  pairs  we  assume   b  >  0 
without  loss  in  generality. 
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Assuming   9   is  complex  allows  us  to  rewrite  (A. 2)  as 


6 (8+e)  =  c(6+e-d) 


k    Mj 


^1=1      1      } 


(A. 3) 


case 


We  consider  the  case   a  ^  0   or   0  <  arg  6  £  y.   In  this 


k    u 
i=l    Ki 


<  1,   and  it  follows  that 


Re 


k   y. 


k   u, 


.*  (eryH-1    <0   or    i<arg Ln  ferpH-1  <f  . 

1=1    1    J  ^1=1   Kl    ; 


Furthermore,   0  <  arg(6+e-d)  <  arg(6+e)  <  arg (6)  £  y,   and 


it  follows  that 


0  <  arg (6 (9+e) )  <  arg 


c(9+e-d) 


ki=l    *i 


<  2lT. 


Consequently,   6   cannot  be  a  root  of  (A. 3)  and  consequently,  (A. 2) 
if   a  ^  0.   We  have  just  proved  that  any  complex  root  must  have 
strictly  negative  real  part,  thus  all   k+1   eigenvalues  of   M 
have  strictly  negative  real  parts. 
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